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Numerical Differentiation: 



In numerical analysis, numerical differentiation describes algorithms for estimating the derivative of 
a mathematical function or function subroutine using values of the function and perhaps other 
knowledge about the function. 

Suppose an Engineer collects a set of data points relating two properties such as temperature and 
pressure, or distance and time that and object has fallen or amount of water that it has received. 

He might wish to determine on or more derivatives of the function but without an analytic 
expression for the function, the derivatives can only be approximated numerically. 



Approximation for the First Derivative: 

Given a function / (x), we can obtain the series expansion of / (x) about x, assuming that the 
function has as many continuous derivatives as may be required . Using Taylor's formula 

/ (x + h) = f (x) + h f\x) + ^ /" (x) + ^ /"' (x) + ... 
Solving for fix) , we get 

f(x + h)-f(x) h h 2 

r(x) = h " 2 r(x) " T r " (x) + ■" 

So that an approximation for fix) may be written as 

fix), 

For small value of h . 

The above expression is known as the Forward Difference Approximation. 

Another approximation for fix) can be obtained by replacing h with -h above, such that 

f(x)-f(x-h) 

/w- - h 

This is known as the Backward Difference Approximation to fix) . 



Note that the truncation error for both the forward and backward difference approximation to 
fix) is — ^ f'ip) , for some p in the interval (x, x + h) , for the forward difference and for 
some p in (x — h , x) , for the backward difference 



If we now place —h instead of h in the Taylor series expansion and then subtract the resulting 
equation from the one above, another form of the approximation to f'(x) is obtained such that 

f(x + K)-f(x-K)=f(x)-f(x) + h /'GO + h /'(*) + ^ f"(.x) - ^ /"(*) + ^ /"'(*) + ^ /"'(*) + •• 
So 

„ N / (x + K)-f (x-h) h 3 „„ N /i 4 ro 
/'GO = J -l J -^-± J - + — /'"(*) - - 

Or 

f{x + h)-f{x-h) 
/( * ) - 2^ 

for a small value of h . 

This is known as the Central Difference Approximation of /'GO , which gives us an approximation 
to /'GO with an error of the order h 2 . 



Approximation for the Second Derivative: 

Considering now approximations to the second derivative /"GO ; consider once more the Taylor 
series expansion of / GO about x 

h 2 h 3 /i (4) 

f{x + h)=f (x) + hf'{x) + — /"GO + - /'"GO + — / W + ... 



And 



/ (x - h ) = / (x) - h f'{x) + ^ /"GO - ^ /'"(x) + ^ /W - ... 



Adding the two equation we get 



2h 2 2h 4 
f(x + h)+f(x-h) = 2f(x) + — /"(*) + 1 ^/ (4) + - 

Now solve for f"(x) to obtain 

f „ ( , /Q + ft)-2/(x)+/(x-/i) 1 

/ GO = ^ ^2 ' ^ 

Therefore, an approximation formula for/"(x) is given by 



For a small value of h . 

In this case it can see that truncation error is — ^ h 2 f^ (p) for some p in the interval (x — 
h,x + h). 

This formula is known as the Central Difference Formula for the second derivative. 



Approximation for the Third Derivative: 

Considering now approximations to the third derivative f"\x), Let's start with Taylor series: 

nx + h) = f(x) + hf\x) + %r\x) + ^/'"(*) + ^j/ (4) + - 

And 

fix -h) = fix) - hf'ix) + - ^rw - ^/ (4) + - 

So the central difference form of first derivative will be 

fix + h) - fix - h) 
/(*)= 2h + '" 

In the previous relation, replace every h with ^ , so the relation will be 



From Taylor series, the second derivative will be 

iff 2h\ f 0h\ f -2h\\ 

The third derivative can be deduced using Taylor expansion: 

f ( x + t) = f w + ! + © 2 i^" w + © 3 S^'" w + - w 

/ (x - y) = /CO - ; V CO + g) 2 y /" CO " g) 3 f f "CO + - 0 

/ (» + s) = m + 1 v'w + © 2 frco + g) 3 £r co + - w 
/ (x - s) = /co -jvw + © 2 Srco - © 3 co + - w 



Subtracting (2) from (1) we get 

f (* + t) - f (*- T) = 2 (f V'M + ©' + ••■ (1 

Subtracting (4) from (3) we get 



f (* + 1) - f {*- 1) - 2 (l vw + ©' I w) + - w 

Multiply (a) by 3 and add it to (b), we get 




It can be written as 



/(-t)-/(-t)-(H-i)^(-i)) = ^'"«- 

So the third derivative will be: 

r w - p (/ (« + 3 4) - v (* 4) + 3 ' (* - s) - ' (* + t)) 

With a truncation error 0(h 2 ) 



Approximation for Higher Order Derivatives: 

It's obvious that from the formula of /', /"and /"' that 

fin) Sf =0 [(-l) 1 * © * / (* + h (| - i))] + r. E. (N) 

Where T.E. is the truncation error always being 0(h 2 ) 
Using this relation we can find that: 



Using Taylor Series get the fourth derivative formula and thus proving the relation (N) 

fix + 2h) = fix) + 2hf\x) + 2 2 |/"(x) + 2 3 ^/'"(x) + 2 4 ^/W + - (5) 

-4/0 + h) = -4 (f(x) + hf'ix) + £/"(*) + g/'"(x) + ^/ (4) + - ) 

-4/(x - h) = -4 (/(x) - fc/'(x) + £/"(*) - g/'"(x) - + ... ) (7j 

/(x - 2/i) = /(x) - 2/i/'(x) + 2 2 |/"(x) - 2 3 |/"'(x) + 2 4 ^/W + - (8) 

By the addition of (5), (6), (7), (8) and 6/(x) we get : 

C.O.fix) = [1-4 + 6-4+1] = 0 
C. 0. hf'ix) = [2-4 + 4-2] = 0 

h 2 

CO. — fix) = [2 2 - 4 - 4 + 2 2 ] = 0 

h 3 

C. 0.—f"'ix) = [2 3 - 4 + 4 - 2 3 ] = 0 
CO. — / (4) (x) = [2 4 - 4-4 + 2 4 ] = 24 

ft 4 

fix + 2h) - 4 fix + h) + 6 fix) - A fix - h) + /(x - 2/i) ^ 24 * -/ (4) (jc) 

4! 

Hence, the fourth derivative is 

/ (4 » M = ^(/(* + 2") " ♦/(* + » + 6/W " Vfr " » + " 2«) + r- «• 

The same idea is applicable to prove (N) for f < - 5 - ) and f ^ and so on. 

/ (5) o) = ^ (f [ x + \ h ) - 5 f [ x + \ h ) + 10 f ( x + \ h ) - 10 f ( x - \ h ) + 5 f [ x - \ h ) - f { x -l h fj 

1 

/ (6) (» = ifix + 3h) - 6 fix + 2h) + 15 fix + h) - 20 fix) + 15 fix - h) - 6 fix - 2h) + fix - 3h)) 



Richardson's Extrapolation: 

The central difference approximation gives a truncation error TE = 0(/i 2 ) , while the forward 
difference approximation gives a truncation error T. E. = 0(h ) . 

We can reduce the T. E. more than 0(/i 2 ) by technique known as "Richardson extrapolation". 

Richardson's extrapolation is an additional technique for approximating derivatives of a function / 
that will enable us to reduce the truncation error. 

We therefore start with the Taylor expansions of f{x ± h) : 

f(x + h)=f(x) + hf f (x) + ^ /"(*) + ^ /"'(*) + 

fix -h)=f(x)- hf'tx) + ^ /"(*) - ^ /"'(*) + ... 
Hence the central difference approximation: 

fix) = f^f^ + £/G» (*) + £ f W (x) + ... (9) 

Rewrite (9) as: 

L = D(h) + c 2 h 2 + c 4 h 4 + .. (10) 
Where L denotes the quantity that we are interested in approximating 

L = f\x) 

And D(h)\s the approximation, which in this case is 

fix + h) - f{x - h) 
D(Jl) = 2h 

The error is 

E = c 2 h 2 + c 4 ft 4 + 
Where q denotes the coefficient of h t in (9). 

By trying to eliminate the term c 2 h 2 the truncation error will be improved. 
If we replace h in equation (10) with ^ 



The terms in h 2 can be eliminated by (4 * (11) - (10)). We get 

3/'(x) = 4D (-) - D(/i) + 4c 2 (-) - c 2 h 2 + 4c 4 (-) - c 4 /i 4 + - 
And it reduces to 

4 /h\ 1 1 
3/'(*) = -Dy- ^D(h)--c 4 h 4 + - 

The accuracy has been improved by having T.E. of 0(/i 4 ) . 
Thus, we get 

4 /h\ 1 
''«-3 D (2)-3 D(A) 

Where / (5 ^ (p) the greatest value of the fifth derivative on interval [x — h , x + h ] and 

/h\ = /(*+f)-/(*-f) . . = f(x+h)-f(x-h) 
\2J h ' y } 2h 

This extrapolation method can be repeated to get higher accuracy by having a 7\ £*. = 0(/i 4 ), 
T.E.= 0(/i 6 )andsoon. 

In general, given and approximation D(K) and having computed the values 

d ul = D{^) f i = l,2,... 
For a given /i > 0, the process can be extended to recursively generate j columns by the formula 

- 47 _ x d U ~ 47 _ x 

The Truncation Error associated with the entry d t j +1 is of order 0(/i 2;+2 ). By arranging the 
quantities in a tabular form, the procedure can be better illustrated 



ij 


h 


d u ,0(/i 2 ) 


di, 2 ,0(h 4 ) 


d ii3 ,0(h 6 ) 




dij +1 ,0(hV +2 ) 


1 


h 


di,i 










2 


h 
2 


d 2 ,i 


d 2 ,2 








3 


h 
4 


d 3 ,i 


d 3 ,2 


d 3 ,3 




















i 


h 


di,i 


d[,2 


di,3 




du 


(2-i) 



Examples: 

1- Consider the function (#) = e x sin(V) . Find the fourth derivative / (4) (2. 2) . 

Solution 

Using h = 0.1 and MATLAB program (first_sixth_derivative.m) 



+-- 

x | 


h 


- + -- 

1 


4 derivative 


+-- 




- + -- 




2.2000 | 


0.1000 


1 


-29.1158 


+ __ 




- + -- 





.-. /W(2.2) = -29.1158 
2- Given the data in the following table: 



t 


0 


0.5 


1.0 


1.5 


2.0 


y 


0 


0.19 


0.26 


0.29 


0.31 



Give the central difference approximations for /"(1),/"'(1),/ (4) (1) . 

Solution 

Using h = 0.5 and MATLAB program (first_sixth_derivative.m) 
- For /"(*) 



+ 





T 






T 






X 


1 


y 




1 


2 


derivative 




+ 






- + 






0 


0000 | 


0 


0000 


I Can 


t 


be calculated 


0 


5000 | 


0 


1900 






-0.4800 


1 


0000 | 


0 


2600 






-0.1600 


1 


5000 | 


0 


2900 






-0.0400 


2 


0000 | 
-j- 


0 


3100 


I Can 


t 


be calculated 



+ 



••/"(I) = -0.1600 



- For /'"(*) 

+ + + + 

I x I Y | 3 derivative | 



+ 





T 






T 








0 


0000 | 


0 


0000 


| Can 


t 


be 


calculated 


0 


5000 | 


0 


1900 


I Can 


t 


be 


calculated 


1 


0000 | 


0 


2600 








0 .4400 


1 


5000 | 


0 


2900 


I Can 


t 


be 


calculated 


2 


0000 | 
-j- 


0 


3100 


I Can 


t 


be 


calculated 



+ 



•/"(l) = 0.4400 

- For/ (4) (x) 

+ + + + 

| x I Y | 4 derivative | 



+ 





T 






T 








0 


0000 | 


0 


0000 


I Can 


t 


be 


calculated 


0 


5000 | 


0 


1900 


I Can 


t 


be 


calculated 


1 


0000 | 


0 


2300 








-1.7600 


1 


5000 | 


0 


2600 


I Can 


t 


be 


calculated 


2 


0000 | 


0 


3100 


I Can 


t 


be 


calculated 




+ 






- + 









+ 



= -1.7600 

- Given the evenly spaced data points, Compute /'(0) 



X 


0 


0.1 


0.2 


0.3 


0.4 


/(*) 


0.0000 


0.0819 


0.1341 


0.1646 


0.1797 



Solution 



Using MATLAB program (first_deriv_forward.m) 



+ + + + 

I x I Y I first derivative | 

+ + + + 

| 0.0000 I 0.0000 I 0.8190 | 
I 0.1000 | 0.0819 | 0.5220 | 
I 0.2000 | 0.1341 | 0.3050 | 
I 0.3000 | 0.1646 | 0.1510 | 
| 0.4000 | 0.1797 | Can s t be calculated | 
+ + + + 



•• /'(0) = 0.8190 



4- In a circuit with impressed voltage e(t) = L— + Ri where: 

at 

R is the resistance in the circuit and i is the current. Suppose we measure the current for 
several values of t and obtain: 



t 


1.00 


1.01 


1.02 


1.03 


1.04 


i 


3.10 


3.12 


3.14 


3.18 


3.24 



Where t is measured in seconds, i is in amperes, the inductance L is a constant 0.98 
henries, and the resistance in 0.142 ohms. Approximate the voltage s(t) when 
t = 1. 00, 1. 01, 1. 02, 1. 03 and . 04 . 

Solution 

Using MATLAB program (first_deriv_forward.m) to get ^ for the first value of t, MATLAB 
program (first_deriv_backward.m) to get ^ for the last value of t, and MATLAB program 
(first_deriv_central.m) to get — for the remaining value of t 



+ + + + 

I x I Y I first derivative | 



+ 





+ 




+ 






1 


0000 | 


3 


1000 | 


2 


0000 


1 


0100 | 


3 


1200 | 


2 


0000 


1 


0200 | 


3 


1400 | 


3 


0000 


1 


0300 | 


3 


1800 | 


5 


0000 


1 


0400 | 
+ 


3 


2400 | 

+ 


6 


0000 



Since = 0.98 H , R = 0.142 ohm 



f(l. 00) = (0. 98)(2. 00) + (0. 142 )(3. 10) = 2. 40702 V (1) 

£(1. 01) = (0. 98)(2. 00) + (0. 142 )(3. 12) = 2. 40304 V (2) 

£(1. 02) = (0. 98)(3. 00) + (0. 142 )(3. 14) = 3. 38588 V (3) 

f(l. 03) = (0. 98)(5. 00) + (0. 142 )(3. 18) = 5. 33156 V (4) 

f(l. 04) = (0. 98)(6. 00) + (0. 142 )(3. 24) = 6. 34008 V (5) 



5- Find approximation to f'(l. 8) for f(x) = \n(x) with h = 0. 05 . Then use Richardson's 
Extrapolation on these values to see if this results in a better approximation. 



Solution 

Using MATLAB program (first_deriv_forward.m) to find /'(1.8) with T.E. = 0(h) 



+ + + + 

I x I h I first derivative | 

+ + + + 

I 1.8000 | 0.0500 | 0.5480 | 
+ + + + 



Using MATLAB program (richardson.m) to use Richardson's Extrapolation with accuracy 
0(h 2 ), 0(/i 4 ), 0(h 6 ) and 0(h s ) 

| i | h | Di,l | Di,2 | Di , 3 | Di,4 | 

I 1 | 0.0500 | 0.5557 I I I I 

I 2 | 0.0250 | 0.5556 | 0.5556 | I I 

I 3 | 0.0125 | 0.5556 | 0.5556 | 0.5556 | I 

I 4 | 0.0063 | 0.5556 | 0.5556 | 0.5556 | 0.5556 | 

Where /'(1.8) = 5555555556 

We see that Richardson's Extrapolation results in a better accuracy lowering the Truncation 
Error. But, it is not very obvious for higher 0(h) because of the Round-Off error. 

6- Taylor Theorem can be used to show that centered-difference formula to approximate 

/'(*o) = ^[/(*o + h) - f(x 0 - h)] - y/'"(* 0 ) " ^/ (5) (*o) - - 

Find approximation with a T.E. of 0(h 2 ), 0(h 4 ) and 0(h 6 ) for /'(2. 0) when f(x) = xe x 
and = 0. 2 . 

Solution 

Using MATLAB program (richardson.m) to use Richardson's Extrapolation with accuracy 
0(h 2 ), 0(h 4 ) and 0(h 6 ) 



I i I h | Di r 1 | Di,2 | Di,3 | 

| 1 | 0.2000 | 22.4142 | | | 

| 2 | 0.1000 | 22.2288 | 22.1670 | | 

| 3 | 0.0500 | 22.1826 | 22.1672 | 22.1672 | 



- Having a T.E. of 0(h 2 ) results in /'(2.0) = 22.4142 

- Having a T.E. of 0(/i 4 ) results in /'(2.0) = 22.1670 

- Having a T.E. of 0(/i 6 ) results in /'(2.0) = 22.1672 

Note that having higher 0(h) results in a better accuracy. 



7- The time t [s] for a body to move from distance 40 meters to a distance d [m] is given 
below for five values of d : 



t[s] 


0.00 


0.50 


1.00 


1.50 


2.00 


d[m] 


40 


43 


60 


77.5 


89 



Find the speed v and acceleration a at t = Is 

Solution 

Since v = — , Using MATLAB program (first_deriv_central.m) to find v (1.00) gives: 



+ + + + 

I x I Y I first derivative | 



+ 





+ 






- + 








0 


.0000 | 


40. 


0000 


I Can s 


t 


be 


calculated 


0 


.5000 | 


43. 


0000 








20.0000 


1 


.0000 | 


60. 


0000 








34.5000 


1 


.5000 | 


77 . 


5000 








29.0000 


2 


.0000 | 
+ 


89. 


0000 


I Can s 
- + 


t 


be 


calculated 



+ 



v (1. 00) = 34. 5 m/s 

Since a = —, Using MATLAB program (first_sixth_derivative.m) to find a (1.00) gives: 

+ + + + 

I x I Y | 2 derivative | 



+ 





+ 






- + 








0 


0000 | 


40 


0000 


| Can 


N t 


be 


calculated 


0 


5000 | 


43 


0000 








56.0000 


1 


0000 | 


60 


0000 








2.0000 


1 


5000 | 


77 


5000 








-24.0000 


2 


0000 | 
+ 


89 


0000 


I Can 
- + 


N t 


be 


calculated 



+ 



.-. a (1. 00) = 2.0 m/s 2 



8- The altitude h [m] attained by an aircraft climbing from sea level at its maximum rate of 
climb [m/s] is shown at four times t [min]: 



t [min] 


10 


20 


30 


40 


h[m] 


5400 


9000 


11600 


12800 



Estimate the aircraft's rate of climb (^) at altitude h = 11600 m 



Solution 



Using MATLAB program (first_deriv_central.m) to calculate /i'(11600)gives 

+ + + + 

I x I Y I first derivative | 

+ + + + 

I 10.0000 I 5400.0000 | Can s t be calculated | 
I 20.0000 | 9000.0000 | 310.0000 | 
I 30.0000 | 11600.0000 | 190.0000 | 
I 40.0000 | 12800.0000 | CarTt be calculated | 
+ + + + 

Hence, aircraft's rate of climb h! at h = 11600 m is 190 m/s. 



9- A rod is rotating in a plane. The following table gives the angle 0 [rad] through which the 
rod has turned for various values of time t [s]. Calculate the angular velocity o) and 
angular acceleration a of the rod at t = 0. 6 s. 



t 


0 


0.2 


0.4 


0.6 


0.8 


1.0 


1.2 


e 


0 


0.12 


0.49 


1.12 


2.02 


3.20 


4.67 



Solution 



d0 

Using MATLAB program (first_sixth_derivative.m) to calculate a) (0.6) = —(0.6) 



1 derivative 





+ 






- + 








0 


0000 | 


0 


0000 


| Can 


N t 


be 


calculated 


0 


2000 | 


0 


1200 








1 .2250 


0 


4000 | 


0 


4900 








2.5000 


0 


6000 | 


1 


1200 








3.8250 


0 


8000 | 


2 


0200 








5.2000 


1 


0000 | 


3 


2000 








6. 6250 


1 


2000 | 
+ 


4 


6700 


| Can 
- + 


N t 


be 


calculated 



co (0.6) = 3. 8250 rad/s 



Using MATLAB program (first_sixth_derivative.m) to calculate a (0.6) = — j(0.6) 



+ + + + 

I x I Y | 2 derivative | 





+ 






- + 








0 


0000 | 


0 


0000 


| Can 


N t 


be 


calculated 


0 


2000 | 


0 


1200 








6.2500 


0 


4000 | 


0 


4900 








6.5000 


0 


6000 | 


1 


1200 








6.7500 


0 


8000 | 


2 


0200 








7.0000 


1 


0000 | 


3 


2000 








7 .2500 


1 


2000 | 
+ 


4 


6700 


| Can 
- + 


N t 


be 


calculated 



/. a (0.6) = 3. 8250 rad/s 2 

10- Given (x) = + tan(e 3x ) . Approximate / (5) (x) and / (6) (x) from x = 0. 5 to 
x = 0.8 with h = 0.02. 



Solution 

Using MATLAB program (n_derivative.m) to for f^(x) gives 



+ + + + 

I x I h I 5 derivative | 
+ + + + 



0 


5000 | 


0 


0200 


34267433518 


8024 


0 


5200 | 


0 


0200 


-42927896685 


2039 


0 


5400 | 


0 


0200 


27380868679 


6787 


0 


5600 | 


0 


0200 


-7845905514 


1852 


0 


5800 | 


0 


0200 


471174460 


9699 


0 


6000 | 


0 


0200 


66660390 


1552 


0 


6200 | 


0 


0200 


172608783 


1919 


0 


6400 | 


0 


0200 


-6345249309 


2422 


0 


6600 | 


0 


0200 


26786279211 


6589 


0 


6800 | 


0 


0200 


-47972667422 


2166 


0 


7000 | 


0 


0200 


43421810751 


8765 


0 


7200 | 


0 


0200 


-19329788330 


8946 


0 


7400 | 


0 


0200 


3608359902 


9256 


0 


7600 | 


0 


0200 


-4064180262 


5284 


0 


7800 | 


0 


0200 


11995601850 


9831 


0 


8000 | 


0 


0200 


-16280002599 


5907 



+ + + + 



Using MATLAB program (n_derivative.m) to for f^(x) gives 

















+ 




+- 






- + 


+ 


1 


X 


1 




h 


| 6 derivative | 
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+- 






- + 


+ 


1 


0 


5000 | 


0 


0200 


| -5785324988021 


5986 | 


1 


0 


5200 | 


0 


0200 


| 7016442104430 


0547 | 


1 


0 


5400 | 


0 


0200 


| -4780347384695 


6123 | 


1 


0 


5600 | 


0 


0200 


| 1700422929929 


8767 | 


1 


0 


5800 | 


0 


0200 


| -230458042287 


5795 | 


1 


0 


6000 | 


0 


0200 


| -3196470565 


1838 | 


1 


0 


6200 | 


0 


0200 


| 42610190635 


8265 | 


1 


0 


6400 | 


0 


0200 


I -423968503056 


7980 | 


1 


0 


6600 | 


0 


0200 


| 1374495883216 


5879 | 




0 


6800 | 


0 


0200 


I -2256680341100 


7251 | 




0 


7000 | 


0 


0200 


| 2072128797821 


9556 | 




0 


7200 | 


0 


0200 


| -1043431143895 


0142 | 




0 


7400 | 


0 


0200 


| -378414453515 


0197 | 




0 


7600 | 


0 


0200 


I 3582703543645 


4082 | 




0 


7800 | 


0 


0200 


I -8732886598788 


3535 | 




0 


8000 | 


0 


0200 


| 11364408063333 


9940 | 
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How MATLAB Programs Work: 

Backward (Forward, Central point) Method: 



Points 
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The input Mode 
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Format 
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Input 




The 




Points 





Input 
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Boundary 



Calculation of 
The 1 st derivatives 

Output of 
The results 



Calculation of 
The 1 st derivatives 



Output of 
The results 



Higher order up to the 6th order: 




Higher Order Derivatives: 




Richardson's Extrapolation: 



Table of 
operations 
for one 
derivative 



Input 
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Points 
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Accuracy 



Choice of 



The output Mode 



Table of many 
derivatives 




Calculation of the 
desired 1st derivative 



Calculation of the 
desired 1st derivatives 



Used MATLAB Functions: 



1- input (parameter!, parameter2): Takes the input from the user 

parameterl : a string parameter. It outputs a massage to tell the user what to input. 
Escape character like "\n", "\\" can be used. 

parameter2 : a character parameter. It is used to when having a string input rather than 
expression. 

2- inline (parameter): Converts string to a function 

parameter : a string parameter that stores the function to a string in order to inline it by the 
inline() function. 

3- disp (parameter): Prints the output on the screen. 
parameter : this can be any type parameter like variable or string. 
The function disp() ends the line after execution. 

You cannot use escape characters like "\n". 

4- strcmp (stringl, string2): Compares two strings and returns logical T if they are identical, 
and logical '0' if they are not. 

stringl & string2: a string parameters. 
The function strcmp() is case insensitive. 

5- abs (parameterl): Returns the absolute value of "paramaterl". 

6- int32 (x): Approximates the value of V to the nearest integer number. 

7- numel (x): Returns the number of elements in V. 

8- fprintf (format, x, y, ...): Format and print the output just like in C++ 

9- double (x): Convert the type of variable V to double. 

10- strcat (stringl, string2): Returns the concatenation of stringl and string2. 
Example: 

stringl = 'TEST1' 
string2 = 'TEST2' 

TEST3 = strcat (stringl, string2) TEST3 = 'TEST1TEST2' 



11- feval (parameter, argl, arg2, arg3, ...): Evaluates the value of "parameter" that has inlined 
content using "parameter" arguments argl, arg2, arg3 ... . 

Example: f = inline ('b/a') -> f(a,b) = b/a 
feval (f, 5, 2) f = 2/5 = 0.4 

12- floor (parameter): Approximates the value of "parameter" to the first lower integer 
number. 



Note: Since MATLAB R2013a warns about removing inline() function in its future release, 
here is another method for inserting a function in MATLAB without using inline() and 
feval() functions. 



x = sym ( 1 x 1 ) ; 


o 
o 


Define a symbol x 


f = input (' Enter f(x)= ! ); 


o 
o 


Enter an "expression" for f (x) 




o 
o 


not "string" as in inline 


f = symfun (f ,x) ; 


o 
o 


Convert sym to sym function 


result = f (5) ; 


o 
o 


Evaluate the function for x = 5 


result = double (result) 


o 
o 


Convert result to double 



First Derivative-Central Difference: 



1. %% Central Difference 

2. % The first derivative " *Central Difference* " metrhod is implemented as 

3. % shown here using MATLAB(R) . 
4. 

5. %% Clear Memory and Screen 

6. clear 

7. clc 
8. 

9. %% Firstly, the user has to choose between two input modes 

10. % 1- Function: Formula of a function. 

11. % 2- Data : some (x,y) points. 

12. % the user has to write one of two words (Function or Data) to indicate the mode. 

13. while 1 

14. disp(' Please Choose your Input mode (Function or Data)'); 

15. Mode = input( ' ' / s ' ) ; 

16. dispC '); 

17. if ~(strcmpi( Mode , 'Function' )||strcmpi( Mode , 'Data' )) 

18. disp( ' Error! ' ) 

19. else 

20. break; 

21. end; 

22. end; 
23. 

24. %% Function Mode 

25. % The user has chosen the Function mode (Go to line 10 ). 

26. if strcmpi( Mode , 'Function' ) 

27. fs = input( ' Please Enter The Function \n f(x)= ','s'); 

28. f = inline(fs); 

29. disp(' '); 

30. 

31. % Here, the user has to define where he wants to get the values of the 1st derivative. 

32. % The user can input that using two formats: points that he wants to get the 1st 
derivative 

33. % or the boundaries(a, b) of an period[a,b] besides a number of periods(n). 

34. while 1 

35. disp(' Please Choose points Input mode (Points or Boundary)'); 

36. dispC Points like [X0 XI X2 ... ]'); 

37. disp( ' Boundary means enter [a,b] and n'); 

38. Point_Mode = input ( ' ' , ' s ' ) ; 

39. disp(' '); 

40. if ~(strcmpi( Point_Mode , 'Boundary' )||strcmpi( Point_Mode , 'Points' )) 

41. disp( ' Error! ' ); 



42. else 

43. break; 

44. end; 

45. end; 

46. % Here, the code deals with the one of the two formats (Go back to Line 31). 

47. if strcmpi( Point_Mode , 'Boundary' ) 

48. a = input( ' Please Enter [ajb]\n a = '); 

49. b = input( ' b = ' ); 

50. n = input( ' Please Enter number of periods\n n = '); 

51. h = (abs(b-a))/n; 

52. X = a:h:b; 

53. h(l:n+l) = h(l); 

54. else 

55. X = input( 'Please Enter X \n For Example X = [0,l,5,6]\n X = '); 

56. % Enter The Value of Step Size h. 

57. disp(' Enter a value of (h) OR the corresponding (h) for every (x)'); 

58. disp('Note: if numbers of (h) is less than number of (x) 3 the remaining'); 

59. disp( ' corresponding values of (h) will be last (h) entered'); 

60. disp( ' '); 

61. h = input( ' Please Enter the value of Step Size\n h = '); 

62. n = numel(X); 

63. nh = numel(h); 

64. % If data not complate we let the other h equal last value of h. 

65. if nh<n 

66. h(nh:n) = h(nh); 

67. end; 

68. end; 

69. disp(' '); 

70. 

71. % Find First Derivative. 

72. n = numel(X); 

73. FD = X; 

74. for i = l:n 

75. FD(i) = (feval(f,X(i)+h(i))-feval(f,X(i)-h(i)))/(2*h(i)); 

76. end; 
77. 

78. %Display X, h, First Derivative. 

79. fprintf ( 'f (x) = %s\n',fs); 

80. disp (' + + + +'); 

81. f printf ( ' | x h | first derivative |\n'); 

82. disp (' + + + +'); 

83. for i = l:n 

84. fprintf (' | %12.4f | %8.4f | %12.4f | \n ' ,X(i) , h(i) , FD(i) ) ; 

85. end; 



86. disp (' + + + +'); 

87. 

88. end; 
89. 

90. %% Data Mode 

91. % The user has chose the Data mode (Go back to line 11 ). 

92. if strcmpi( Mode , 'Data' ) 



93. % The user has to enter the Boundary of interval. 

94. a = input( ' Please Enter the interval [a,b]\n a = '); 

95. b = inputC b = '); 

96. n = input( ' Please Enter number of periods\n n = ' ); 

97. if b<a 

98. tmp = a; 

99. a = b; 

100. b = tmp; 

101. end; 

102. h = (b-a)/n; 

103. X = a:h:b; 
104. 

105. % The user has to enter Y(x). 

106. Y = X; 

107. for i = l:n+l 

108. fprintf( 'Please Enter Y(%1.4f) = 1 ,X(i)); 

109. Y(i) = input( 1 ' ); 

110. end; 
111. 

112. % Find First Derivative. 

113. FD = X; 

114. for i = 2:n 

115. FD(i) = (Y(i+l)-Y(i-l))/(2*h); 

116. end; 
117. 

118. %Display X, Y, First Derivative. 

119. disp (' + + + +'); 

120. f printf ( ' | x y first derivative |\n'); 

121. disp (' + + + +'); 

122. fprintf('| %12.4f | %12.4f | Can^t be calculated | \n ' ,X(1),Y(1)); 

123. for i = 2:n 

124. fprintf('| %12.4f | %12.4f | %12.4f | \n ' ,X(i) , Y(i) , FD(i) ) ; 

125. end; 

126. fprintf('| %12.4f | %12.4f | Can^t be calculated | \n ' ,X(n+l),Y(n+l)); 

127. disp (' + + + +'); 

128. 

129. end; 



First Derivative-Forward Difference: 



1. %% Forward Difference 

2. % The first derivative " *Forward Difference* " metrhod is implemented as 

3. % shown here using MATLAB(R) . 
4. 

5. %% Clear Memory and Screen 

6. clear 

7. clc 
8. 

9. %% Firstly, the user has to choose between two input modes 

10. % 1- Function: Formula of a function. 

11. % 2- Data : some (x,y) points. 

12. % the user has to write one of two words (Function or Data) to indicate the mode. 

13. while 1 

14. disp(' Please Choose your Input mode (Function or Data)'); 

15. Mode = input( ' ' / s ' ) ; 

16. dispC '); 

17. if ~(strcmpi( Mode , 'Function' )||strcmpi( Mode , 'Data' )) 

18. disp( ' Error! ' ) 

19. else 

20. break; 

21. end; 

22. end; 
23. 

24. %% Function Mode 

25. % The user has chosen the Function mode (Go to line 10 ). 

26. if strcmpi( Mode , 'Function' ) 

27. fs = input( ' Please Enter The Function \n f(x)= ','s'); 

28. f = inline(fs); 

29. disp(' '); 

30. 

31. % Here, the user has to define where he wants to get the values of the 1st derivative. 

32. % The user can input that using two formats: points that he wants to get the 1st 
derivative 

33. % or the boundaries(a, b) of an period[a,b] besides a number of periods(n). 

34. while 1 

35. disp(' Please Choose points Input mode (Points or Boundary)'); 

36. disp(' Points like [X0 XI X2 ... ]'); 

37. disp( ' Boundary means enter [a,b] and n'); 

38. Point_Mode = input ( ' ' , ' s ' ) ; 

39. disp(' '); 

40. if ~(strcmpi( Point_Mode , 'Boundary' )||strcmpi( Point_Mode , 'Points' )) 

41. disp( ' Error! ' ); 



42. else 

43. break; 

44. end; 

45. end; 

46. % Here, the code deals with the one of the two formats (Go back to Line 31). 

47. if strcmpi( Point_Mode , 'Boundary' ) 

48. a = input( ' Please Enter [ajb]\n a = '); 

49. b = input( ' b = ' ); 

50. n = input( ' Please Enter number of periods\n n = '); 

51. h = (abs(b-a))/n; 

52. X = a:h:b; 

53. h(l:n+l) = h(l); 

54. else 

55. X = input( 'Please Enter X \n For Example X = [0,l,5,6]\n X = '); 

56. % Enter The Value of Step Size h. 

57. disp(' Enter a value of (h) OR the corresponding (h) for every (x)'); 

58. disp('Note: if numbers of (h) is less than number of (x) 3 the remaining'); 

59. disp( ' corresponding values of (h) will be last (h) entered'); 

60. disp( ' '); 

61. h = input( ' Please Enter the value of Step Size\n h = '); 

62. n = numel(X); 

63. nh = numel(h); 

64. % If data not complate we let the other h equal last value of h. 

65. if nh<n 

66. h(nh:n) = h(nh); 

67. end; 

68. end; 

69. disp(' '); 

70. 

71. % Find First Derivative. 

72. n = numel(X); 

73. FD = X; 

74. for i = l:n 

75. FD(i) = (feval(f,X(i)+h(i))-feval(f J X(i)))/(h(i)); 

76. end; 
77. 

78. %Display X, h, First Derivative. 

79. fprintf ( 'f (x) = %s\n',fs); 

80. disp (' + + + +'); 

81. f printf ( ' | x h | first derivative |\n'); 

82. disp (' + + + +'); 

83. for i = l:n 

84. fprintf (' | %12.4f | %8.4f | %12.4f | \n ' ,X(i) , h(i) , FD(i) ) ; 

85. end; 



86. 


disp (' + + + + '); 


87. 






88 . 


end; 




89. 






90. 


%% Data 


Mode 


01 

yi . 


% The user has chose the Data mode (Go back to line 11 ). 


m 

yz . 


if strcmpi( Mode , 'Data' ) 


93 . 


% 


The user has to enter the Boundary of interval. 


y4 . 


a 


= input( ' Please Enter the interval [a,b]\n a = '); 


95 . 


b 


= input( ' b = ' ) ; 


y© . 


n 


= input( ' Please Enter number of periods\n n = ' ); 


97 . 


if 


b<a 


y© . 




tmp = a; 


99 . 




a = b; 


100. 




b = tmp; 


101 . 




end; 


102. 




h = (b-a)/n; 


103. 




X = a:h:b; 


104. 






1 f>C 

10b . 




% The user has to enter Y(x). 


106. 




Y = X; 


107. 




for i = l:n+l 


108. 




fprintf( 'Please Enter Y(%1.4f) = 1 ,X(i)); 


109 . 




Y(i) = input("); 


110. 




end; 


111 . 






112 . 




% Find First Derivative. 


113 . 




FD = X; 


114. 




for i = l:n 


115 . 




FD(i) = (Y(i + l)-Y(i))/h; 


116. 




end; 


1 1 "7 
11/ . 






118. 




%Display X, Y, First Derivative. 


119 . 




disp (' + + + +'); 


120. 




f printf ( ' | x y first derivative |\n'); 


<i T| 

1Z1 . 




disp (' + + + +'); 


IZz . 




for i = l:n 


123 . 




fprintfC | %12.4f |%12.4f | %12.4f | \n ' ,X(i) , Y(i) , FD(i) ) ; 


124. 




end; 


125. 




fprintfC | %12.4f | %12.4f | Can^t be calculated | \n ' ,X(n+l) , Y(n+1) ) ; 


126. 




disp (' + + + +'); 


127. 






128. 






129. 




end; 



First Derivative-Backward Difference: 



1. %% Backward Difference 

2. % The first derivative " *% Backward Difference Difference* " metrhod is implemented as 

3. % shown here using MATLAB(R) . 
4. 

5. %% Clear Memory and Screen 

6. clear 

7. clc 
8. 

9. %% Firstly, the user has to choose between two input modes 

10. % 1- Function: Formula of a function. 

11. % 2- Data : some (x,y) points. 

12. % the user has to write one of two words (Function or Data) to indicate the mode. 

13. while 1 

14. disp(' Please Choose your Input mode (Function or Data)'); 

15. Mode = input( ' ' / s ' ) ; 

16. dispC '); 

17. if ~(strcmpi( Mode , 'Function' )||strcmpi( Mode , 'Data' )) 

18. disp( ' Error! ' ) 

19. else 

20. break; 

21. end; 

22. end; 
23. 

24. %% Function Mode 

25. % The user has chosen the Function mode (Go to line 10 ). 

26. if strcmpi( Mode , 'Function' ) 

27. fs = input( ' Please Enter The Function \n f(x)= ','s'); 

28. f = inline(fs); 

29. disp(' '); 

30. 

31. % Here, the user has to define where he wants to get the values of the 1st derivative. 

32. % The user can input that using two formats: points that he wants to get the 1st 
derivative 

33. % or the boundaries(a, b) of an period[a,b] besides a number of periods(n). 

34. while 1 

35. disp(' Please Choose points Input mode (Points or Boundary)'); 

36. dispC Points like [X0 XI X2 ... ]'); 

37. disp( ' Boundary means enter [a,b] and n'); 

38. Point_Mode = input ( ' ' , ' s ' ) ; 

39. disp(' '); 

40. if ~(strcmpi( Point_Mode , 'Boundary' )||strcmpi( Point_Mode , 'Points' )) 

41. disp( ' Error! ' ); 



42. else 

43. break; 

44. end; 

45. end; 

46. % Here, the code deals with the one of the two formats (Go back to Line 31). 

47. if strcmpi( Point_Mode , 'Boundary' ) 

48. a = input( ' Please Enter [ajb]\n a = '); 

49. b = input( ' b = ' ); 

50. n = input( ' Please Enter number of periods\n n = '); 

51. h = (abs(b-a))/n; 

52. X = a:h:b; 

53. h(l:n+l) = h(l); 

54. else 

55. X = input( 'Please Enter X \n For Example X = [0,l,5,6]\n X = '); 

56. % Enter The Value of Step Size h. 

57. disp(' Enter a value of (h) OR the corresponding (h) for every (x)'); 

58. disp('Note: if numbers of (h) is less than number of (x) 3 the remaining'); 

59. disp( ' corresponding values of (h) will be last (h) entered'); 

60. disp( ' '); 

61. h = input( ' Please Enter the value of Step Size\n h = '); 

62. n = numel(X); 

63. nh = numel(h); 

64. % If data not complate we let the other h equal last value of h. 

65. if nh<n 

66. h(nh:n) = h(nh); 

67. end; 

68. end; 

69. disp(' '); 

70. 

71. % Find First Derivative. 

72. n = numel(X); 

73. FD = X; 

74. for i = l:n 

75. FD(i) = (feval(f,X(i))-feval(f,X(i)-h(i)))/(h(i)); 

76. end; 
77. 

78. %Display X, h, First Derivative. 

79. fprintf ( 'f (x) = %s\n',fs); 

80. disp (' + + + +'); 

81. f printf ( ' | x h | first derivative |\n'); 

82. disp (' + + + +'); 

83. for i = l:n 

84. fprintf (' | %12.4f | %8.4f | %12.4f | \n ' ,X(i) , h(i) , FD(i) ) ; 

85. end; 



86. disp (' + + + +'); 

87. 

88. end; 
89. 

90. %% Data Mode 

91. % The user has chose the Data mode (Go back to line 11 ). 

92. if strcmpi( Mode , 'Data' ) 



93. % The user has to enter the Boundary of interval. 

94. a = input( ' Please Enter the interval [a,b]\n a = '); 

95. b = inputC b = '); 

96. n = input( ' Please Enter number of periods\n n = ' ); 

97. if b<a 

98. tmp = a; 

99. a = b; 

100. b = tmp; 

101. end; 

102. h = (b-a)/n; 

103. X = a:h:b; 
104. 

105. % The user has to enter Y(x). 

106. Y = X; 

107. for i = l:n+l 

108. fprintf( 'Please Enter Y(%1.4f) = 1 ,X(i)); 

109. Y(i) = input( 1 ' ); 

110. end; 
111. 

112. % Find First Derivative. 

113. FD = X; 

114. for i = 2:n+l 

115. FD(i) = (Y(i)-Y(i-l))/(h); 

116. end; 
117. 

118. %Display X, Y, First Derivative. 

119. disp (' + + + +'); 

120. f printf ( ' | x y first derivative |\n'); 

121. disp (' + + + +'); 

122. fprintf('| %12.4f | %12.4f | Can^t be calculated | \n ' ,X(1),Y(1)); 

123. for i = 2:n+l 

124. fprintf('| %12.4f | %12.4f | %12.4f | \n ' ,X(i) , Y(i) , FD(i) ) ; 

125. end; 

126. disp (' + + + +'); 

127. 

128. 

129. end; 



First To Sixth Order Derivative-Central Difference: 



1. % Central Difference 

2. % The Higher derivative Up to 6 " *Central Difference* " metrhod is implemented as 

3. % shown here using MATLAB(R) . 
4. 

5. %% Clear Memory and Screen 

6. clear 

7. clc 
8. 

9. %% Firstly, the user has to choose between two input modes 

10. % 1- Function: Formula of a function. 

11. % 2- Data : some (x,y) points. 

12. % the user has to write one of two words (Function or Data) to indicate the mode. 

13. while 1 

14. disp(' Please Choose your Input mode (Function or Data)'); 

15. Mode = input( ' ' / s ' ) ; 

16. dispC '); 

17. 

18. if ~(strcmpi( Mode , 'Function' )||strcmpi( Mode , 'Data' )) 

19. disp( ' Error! ' ) 

20. else 

21. break; 

22. end; 

23. end; 

24. %% Enter the order of Derivative 

25. while 1 

26. disp(' Please Enter The order of Derivative (1 to 6)'); 

27. Order = input(' Order = '); 

28. disp(' '); 

29. if (Order ~= l)&&(Order ~= 2)&&(0rder ~= 3)&&(0rder ~= 4)&&(0rder ~= 5)&&(0rder ~= 6) 

30. disp( ' Error! ' ); 

31. else 

32. break; 

33. end; 

34. end; 
35. 

36. %% Function Mode 

37. % The user has chosen the Function mode (Go to line 10 ). 

38. if strcmpi( Mode , 'Function' ) 

39. % Enter the function 

40. fs = input( ' Please Enter The Function \n f(x)= ','s'); 

41. f = inline(fs); 

42. disp(' '); 



43. 

44. % Here, the user has to define where he wants to get the values of the 1st derivative. 

45. % The user can input that using two formats: points that he wants to get the 1st 
derivative 

46. % or the boundaries(a s b) of an period[ajb] besides a number of periods(n). 

47. while 1 

48. disp(" Please Choose points Input mode (Points or Boundary)'); 

49. dispC Points like [X0 XI X2 ... ]'); 

50. disp( ' Boundary means enter [a,b] and n'); 

51. Point_Mode = input ( ' ' , ' s ' ) ; 

52. dispC '); 

53. if ~(strcmpi( Point_Mode , 'Boundary' )||strcmpi( Point_Mode , 'Points' )) 

54. disp( ' Error! ' ); 

55. else 

56. break; 

57. end; 

58. end; 
59. 

60. % Here, the code deals with the one of the two formats (Go back to Line 44). 

61. if strcmpi( Point_Mode , 'Boundary' ) 

62. a = input( ' Please Enter [ajb]\n a = '); 

63. b = input( ' b = ' ); 

64. n = input( ' Please Enter number of periods\n n = '); 

65. h = (abs(b-a))/n; 

66. X = a:h:b; 

67. h(l:n+l) = h(l); 

68. else 

69. X = input( 'Please Enter X \n For Example X = [0,l,5,6]\n X = '); 

70. %Enter The Value of Step Size h 

71. disp(' Enter a value of (h) OR the corresponding (h) for every (x)'); 

72. disp('Note: if numbers of (h) is less than number of (x), the remaining'); 

73. disp( ' corresponding values of (h) will be last (h) entered'); 

74. disp c ■). 

75. h = input( ' Please Enter the value of Step Size\n h = '); 

76. n = numel(X); 

77. nh = numel(h); 

78. % If data not complate for unequal we let the other h equal last value of h 

79. if nh<n 

80. h(nh:n) = h(nh); 

81. end; 

82. end; 

83. disp(' '); 

84. 

85. % Find Needed Derivative. 



86. n = numel(X); 

87. HD = X; 

88. for i = l:n 

89. if Order == 1 

90. HD(i) = (feval(f J X(i)+h(i))-feval(f ,X(i)-h(i)))/(2*h(i)); 

91. elseif Order == 2 

92. HD(i) = (feval(f,X(i)+h(i))+feval(f,X(i)-h(i))-2*feval(f,X(i)))/(h(i)*h(i)); 

93. elseif Order == 3 

94. HD(i) = (feval(f,X(i)+2*h(i))-2*feval(f,X(i)+h(i))+2*feval(f,X(i)-h(i))-feval(f,X(i)- 
2*h(i)))/(2*h(i) A 3); 

95. elseif Order == 4 

96. HD(i) = (feval(f,X(i)+2*h(i))-4*feval(f,X(i)+h(i))+6*feval(f,X(i))-4*feval(f,X(i)- 
h(i))+feval(f,X(i)-2*h(i)))/(h(i) A 4); 

97. elseif Order == 5 

98. HD(i) = (feval(f,X(i)+3*h(i))-4*feval(f,X(i)+2*h(i))+5*feval(f,X(i)+h(i))- 
5*feval(f J X(i)-h(i))+4*feval(f,X(i)-2*h(i))-feval(f,X(i)-3*h(i)))/(2*h(i) A 5); 

99. elseif Order == 6 

100. HD(i) = (feval(f,X(i)+3*h(i))-6*feval(f,X(i)+2*h(i))+15*feval(f,X(i)+h(i))- 
20*feval(f ,X(i))+15*feval(f,X(i)-h(i))-6*feval(f ,X(i) -2*h(i) )+feval(f ,X(i) -3*h(i) ) )/(h(i) A 6) ; 

101. end; 

102. end; 
103. 

104. %Display X, Needed Derivative. 

105. fprintf ( 'f (x) = %s\n',fs); 

106. disp (' + + + +'); 

107. fprintf('| x | h | %1.0f derivative |\n', Order); 

108. disp (' + + + +'); 

109. for i = l:n 

110. fprintf (' | %12.4f | %8.4f | %12.4f | \n ' ,X(i) , h(i) , HD(i) ) ; 

111. end; 

112. disp (' + + + +'); 

113. end; 
114. 

115. %% Data Mode 

116. % The user has chose the Data mode (Go back to line 11 ). 

117. if strcmpi( Mode y 'Data' ) 

118. % The user has to enter the Boundary of interval. 

119. a = input( ' Please Enter the interval [a^bJXn a = '); 

120. b = input(" b = '); 

121. n = input( ' Please Enter number of periods\n n = '); 

122. if b<a 

123. tmp = a; 

124. a = b; 

125. b = tmp; 



126. end; 

127. h = (b-a)/n; 

128. X = a:h:b; 
129. 

130. % The user has to enter Y(x). 

131. Y = X; 

132. for i = l:n+l 

133. fprintf( 'Please Enter Y(%1.4f) = ',X(i)); 

134. Y(i) = input(' ' ); 

135. end; 
136. 

137. % Find Needed Derivative. 

138. n = numel(X); 

139. HD = X; 

140. HP = floor((0rder+l)/2); % half of points that we need to find the drivative 

141. for i = HP+1 : n-HP 

142. if Order == 1 

143. HD(i) = (Y(i+l)-Y(i-l))/(2*h); 

144. elseif Order == 2 

145. HD(i) = (Y(i+l)+Y(i-l)-2*Y(i))/(h*h); 

146. elseif Order == 3 

147. HD(i) = (Y(i+2)-2*Y(i+l)+2*Y(i-l)-Y(i-2))/(2*h A 3); 

148. elseif Order == 4 

149. HD(i) = (Y(i+2)-4*Y(i+l)+6*Y(i)-4*Y(i-l)+Y(i-2))/(h A 4); 

150. elseif Order == 5 

151. HD(i) = (Y(i+3)-4*Y(i+2)+5*Y(i+l)-5*Y(i-l)+4*Y(i-2)-Y(i-3))/(2*h A 5); 

152. elseif Order == 6 

153. HD(i) = (Y(i+3)-6*Y(i+2)+15*Y(i+l)-20*Y(i)+15*Y(i-l)-6*Y(i-2)+Y(i-3))/(h A 6); 

154. end; 

155. end; 
156. 

157. %Display X, Y, H, HD 

158. disp (' + + + +'); 

159. fprintf('| x | y | %1.0f derivative |\n', Order); 

160. disp (' + + + +'); 

161. for i = 1 : n 

162. if (i<HP+l) | | (i>n-HP) 

163. fprintf('| %12.4f | %12.4f | Can^t be calculated | \n ' ,X(i) ,Y(i) ) ; 

164. else 

165. fprintf('| %12.4f | %12.4f | %12.4f | \n ' ,X(i),Y(i),HD(i)); 

166. end; 

167. end; 

168. disp (' + + + +'); 



Higher Order Derivative-Central Difference: 



1. %% Central Difference 

2. % The Higher derivative Up to n " *Central Difference* " metrhod is implemented as 

3. % shown here using MATLAB(R) . 
4. 

5. %% Clear Memory and Screen 

6. clear 

7. clc 
8. 

9. %% Here, The user enter the order of derivative that he needs. 

10. disp(' Please Enter The order of Derivative'); 

11. Order = input(' Order = '); 

12. disp(' ■); 

13. 

14. %% Enter The Function. 

15. fs = input( ' Please Enter The Function \n f(x)= ','s'); 

16. f = inline(fs); 
17. 

18. %% Input Mode 

19. % Here, the user has to define where he wants to get the values of the 1st derivative. 

20. % The user can input that using two formats: points that he wants to get the 1st derivative 

21. % or the boundaries^, b) of an period[a,b] besides a number of periods(n). 

22. disp(' '); 



23. while 1 



29. 



24. 



25. 



26. 



27. 



28. 



disp(' Please Choose points Input mode (Points or Boundary)'); 
disp(' Points like [X0 XI X2 ... ]'); 
disp( ' Boundary means enter [a,b] and n'); 
Point_Mode = input( ' ' , ' s ' ) ; 

disp(' '); 

if ~(strcmpi( Point_Mode , 'Boundary' )||strcmpi( Point_Mode , 'Points' )) 



30. 



disp( ' Error! ' ); 



31. 



else 



32. 



break; 



33. 



end; 



34. end; 

35. % Here, the code deals with the one of the two formats (Go back to Line 19). 

36. if strcmpi( Point_Mode , 'Boundary' ) 

37. a = input( ' Please Enter [a,b]\n a = '); 

38. b = input( ' b = ' ); 

39. n = input( ' Please Enter number of periods\n n = '); 

40. h = (abs(b-a))/n; 

41. X = a:h:b; 

42. h(l:n+l) = h(l); 



43. else 

44. X = input( 'Please Enter X \n For Example X = [0,1,5,61X11 X = '); 

45. %Enter The Value of Step Size h. 

46. disp(' Enter a value of (h) OR the corresponding (h) for every (x)'); 

47. disp('Note: if numbers of (h) is less than number of (x), the remaining'); 

48. disp( ' corresponding values of (h) will be last (h) entered'); 

49. disp(' '); 

50. h = input( ' Please Enter the value of Step Size\n h = '); 

51. n = numel(X); 

52. nh = numel(h); 

53. % If data not complate for unequal we let the other h equal last value of h. 

54. if nh<n 

55. h(nh:n) = h(nh); 

56. end; 

57. end; 

58. disp(' '); 

59. 

60. % Find the Needed Derivative. 

61. n = numel(X); 

62. HD = X; 

63. for j = l:n 

64. HD(j) = 0; 

65. for i = 0: Order 

66. HD ( j ) = HD ( j ) + ((-l) A i) * nchoosek(Order, i) * feval(f ,X(j )+(( (Order/2) -i)*h(j ))) ; 

67. end; 

68. HD( j ) = HD(j)/(h(j) A Order); 

69. end; 
70. 

71. %Display X, h, HD 

72. fprintf ( 'f (x) = %s\n',fs); 

73. disp ( '+ + + +' ); 

74. fprintf (' | x | h | %1.0f derivative | \n ' , Order); 

75. disp ( '+ + + +' ); 

76. for i = l:n 

77. fprintf (' | %12.4f | %8.4f | %12.4f | \n ' ,X(i) , h(i) , HD(i) ) ; 

78. end; 

79. disp ( '+ + + +' ); 



Richardson's Extrapolation: 



1. %% Richardson Extrapolation 

2. % The first derivative " *Richardson Extrapolation* " method is implemented as 

3. % shown here using MATLAB(R) . 
4. 

5. %% Clear Momory and Screen 

6. clear 

7. clc 

8. %% Enter The Function 

9. fs = input( ' Please Enter The Function \n f(x)= ','s'); 

10. f = inline(fs); 

11. dispC '); 

12. 

13. %% Choose Output Mode 

14. % 1- Direvative table (Dnl, Dn2, Dn3, ) for only one point. 

15. % 2- only the value of Derivative for many value of x. 

16. while 1 

17. disp(' choose output mode'); 

18. disp('l => Direvative table (Dnl, Dn2, Dn3, ....) for only one point'); 

19. disp('2 => only the value of Derivative for many value of x'); 

20. Out_Mode = input(' Enter mode (1 or 2) = '); 

21. dispC '); 

22. if Out_Mode == 1 | | Out_Mode == 2 

23. break; 

24. else 

25. disp( ' Error ' ); 

26. end; 

27. end; 
28. 

29. %% Input Mode 

30. % Here, the user has to define where he wants to get the values of the 1st derivative. 

31. % The user can input that using two formats: points that he wants to get the 1st derivative 

32. % or the boundaries^, b) of an period[a,b] besides a number of periods(n). 
33. 

34. % 1- Direvative Table Mode. 

35. if Out_Mode == 1 

36. Point_Mode = 'Points'; %Dust one Point X(l). 

37. end; 

38. % 2- Value Mode. 

39. while Out_Mode == 2 

40. disp(' Please Choose points Input mode (Points or Boundary)'); 

41. disp(" Points like [X0 XI X2 ... ]'); 

42. disp( ' Boundary means enter [a,b] and n'); 



43. Point_Mode = input ( ' ' y ' s ' ) ; 

44. disp c - ); 

45. if ~(strcmpi( Point_Mode , 'Boundary' )||strcmpi( Point_Mode , 'Points' )) 

46. disp( ' Error! ' ); 

47. else 

48. break; 

49. end; 

50. end; 
51. 

52. %% Enter Boundary or Points :: Values of X 

53. if strcmpi( Point_Mode , 'Boundary' ) 



54. a = input( ' Please Enter [a,b]\n a = '); 

55. b = inputC b = '); 

56. n = input( ' Please Enter number of periods\n n = '); 

57. h = (abs(b-a))/n; 

58. X = a:h:b; 

59. h(l:n+l) = h(l); 

60. else 

61. X = input( ' Please Enter the value of x = '); 

62. end; 

63. disp(' '); 

64. 



65. %% Choose the accuracy (h^m) 

66. disp( 'Choose the accuracy you need'); 

67. disp('If the T.E = 0(h A m) and h<l then for greater m we get small error'); 

68. while 1 

69. disp('Note: m must be even number'); 

70. m = input( 'Enter m = '); 

71. disp(' '); 

72. if (mod(m,2) == 0) 

73. break; 

74. end; 

75. end; 

76. m = int32(m/2); 
77. 

78. %% Enter The Value of Step Size (h) for points mode 

79. if strcmpi( Point_Mode , 'Points' ) 

80. if Out_Mode == 2 

81. disp(' Enter a value of (h) OR the corresponding (h) for every (x)'); 

82. disp('Note: if numbers of (h) is less than number of (x), the remaining'); 

83. disp( ' corresponding values of (h) will be last (h) entered'); 

84. disp(' '); 

85. end; 

86. h = input( ' Please Enter the value of Step Size\n h = '); 



87. 


disp(' '); 


oo 
oo . 


n = 


numel(X); 


89. 


nh 


= numel(h); 


90 . 


% If data not complate for unequal we let the other h equal last value of h. 


91 . 


if 


nh<n 


92 . 




h(nh:n) = h(nh); 


93 . 


end; 


y4. end 


y 




95. 






96. %% 


Find 


First Drivative 


97. 






98. D = 


zeros(m); % prepare size of the matrix D to avoid resizing. 


99. % First 


Derivative of X(j) ==> FD 


100. 




FD = X; % prepare size of the matrix FD to avoid resizing. 


101. 






102. 




% Direvative Table Mode just use X(l). 


103. 




% This Condition to deal with more than one value of X. 


104. 




if Out_Mode == 1 


105. 




n = 1; 


106. 




else 


107. 




n = numel(X); 


108. 




end; 


109. 






110. 




for j = 1 : n 


111. 




hj = h(j); 


112. 




D(l,l) = (feval(f,X(j) + hj)-feval(f,X(j)-hj))/(2*hj); 


113. 




for i = l:l:m-l 


114. 




hj = hj/2; 


115. 




D(i+l,l) = (feval(f,X(j)+hj)-feval(f J X(j)-hj))/(2*hj); 


116. 




for k = 1:1: i 


117. 




D(i+l,k+l) = D(i+l,k)+(D(i+l,k)-D(i,k))/((4.0 A double(k))-l); 


118. 




end; 


119. 




end; 


120. 




FD(j) = D(m,m); 


121. 




end; 


122. 






123. 




%% Output for Table Mode 


124. 




if Out_Mode == 1 


125. 




fprintf('f(x) = %s\n',fs); 


126. 




Line = '+ + +' ; 


127. 




for i = l:m 


128. 




Line = strcat ( Line, ' +'); 


129. 




end; 


130. 







131. 


disp(Line); % Before Header. 






liz . 


f printf ( | i | h | ), 






133 . 


for i = l:m 






134. 


tprintt( Dij7ol.0d 


1 ' A \ • 




135. 


end; 






136 . 


f printf ( ' \n ' ) ; 






137. 








1 O o 


disp(Line); % Before Data. 






1 DO 


for i = l:m 






140. 


r~ • i_f"/i| o/ /~v c o/(-i /ir" 

fprintf( | %2.0f %8.4f 


'jijh); 




1/11 

141 . 


f o r j = 1 : m 






142 . 


if j>i 






143 . 


f printf ( 1 


i ' \ . 

1 ); 




144. 


else 






145 . 


tprintt( 7olz.4t | 


r\ / A -A \ \ . 




146 . 


end; 






147. 


end; 






148 . 


f printf ( ' \n ' ) ; 






149 . 


n =h/2; 






150. 


end; 






151 . 








152. 


disp(Line); % End. 






1 C "D 

lbo . 


end ; 






1 C A 

1d4 . 








155 . 


/o/o Output for Value Mode 






156 . 


if Out_Mode == 2 






157. 


fprintf( f(x) = %s\n ,fs); 






158 . 


disp ( ' + + 


1 


+ y jj 


159. 


fprintf( | x h 


| first 


derivative |\n'); 


1 CO 

150 . 


disp ( + + 


1 


+ )> 


161. 


for i = l:n 






162. 


fprintfC | %12.4f | %8.4f 


%12.4f 


|\n , J X(i) J h(i),FD(i)); 


163. 


end; 






164. 


disp (' + + 


+ 


+'); 


165. 


end; 
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